Multiplet effects in orbital and spin ordering phenomena: A hybridization-expansion 

quantum impurity solver study 



Andreas Flesch and Evgeny Gorelov 
Institute for Advanced Simulation, Forschungszentrum Jiilich, 52425 Jiilich, Germany 



Erik Koch 

German Research School for Simulation Science, Jiilich, 
JARA High-Performance Computing 



Germany and 



Eva Pavarini 

Institute for Advanced Simulation, Forschungszentrum Jiilich, 52425 Jiilich, Germany and 

JARA High-Performance Computing 

Orbital and spin ordering phenomena in strongly correlated systems are studied using the local- 
density approximation + dynamical mean-field theory approach. Typically, however, such simu- 
lations are restricted to simplified models (density-density Coulomb interactions, high symmetry 
couplings and few-band models). In this work we implement an efficient general hybridization- 
expansion continuous-time quantum Monte Carlo impurity solver (Krylov approach) which allows 
us to investigate orbital and spin ordering in a more realistic setting, including interactions that 
are often neglected (e.g., spin-flip and pair- hopping terms), enlarged basis sets (full d versus e g ), 
low-symmetry distortions, and reaching the very low-temperature (experimental) regime. We use 
this solver to study ordering phenomena in a selection of exemplary low-symmetry transition-metal 
oxides: LaMn03 and rare-earth manganites as well as the perovskites CaV03 and YTi03. We find 
that, in all considered cases, the minus sign problem mostly appears when off-diagonal crystal-field 
terms are present (and is strongly suppressed in the basis of crystal-field states), while off-diagonal 
terms of the hybridization function matrix are not as critical. We show that spin-flip and pair 
hopping terms do not affect the Kugel-Khomskii orbital-order melting transition in rare-earth man- 
ganites, or the suppression of orbital fluctuations driven by crystal field and Coulomb repulsion. For 
the Mott insulator YTiOs we find a ferromagnetic transition temperature Tc ~ 50 K, in remarkably 
good agreement with experiments. For LaMnOs we show that the classical f2 S -spin approxima- 
tion, commonly adopted for studying manganites, yields indeed an occupied e g orbital in very good 
agreement with that obtained for the full d 5-orbital Hubbard model, while the spin-spin e g -t2 g 
correlation function calculated from the full d model is ~ 0.74, very close to the value expected for 
aligned e g and ti g spins; the e g spectral function matrix is also well reproduced. Finally, we show 
that the t2 g screening reduces the e g -e g Coulomb repulsion by about 10%. 

PACS numbers: 71.10.Fd, 71.10.-w,71.27.+a,71.28.+d,71.30.+h 



I. INTRODUCTION 



Orbital and magnetic ordering phenomena play a cru- 
cial role in the physics of strongly correlated transition- 
metal oxides. Their onset depends on symmetry, lat- 
tice distortions, super-exchange interaction and the form 
of the Coulomb tensor. The realistic description of or- 
dering phenomena requires the ability of disentangling 
the effects of all these interactions. In recent years, 
the local-density approximation+dynamical mean-field 
theory approach^— (LDA+DMFT), which combines ab- 
initio techniques based on density functional theory in 
the local-density approximation (LDA), and the dynam- 
ical mean- field theory^ (DMFT), has lead to important 
progress in understanding such ordering phenomena. It 
has been shown that many-body super-exchange only 
weakly affects the onset of the orbital-order to disorder 
transition in rare-earth manganites,— while, in the pres- 
ence of strong Coulomb repulsion, a small crystal-field 
is sufficient to strongly suppress orbital fluctuations and 



stabilize orbital order However, the effects of subtle 
Coulomb interactions, such as spin-flip and pair-hopping 
terms or of quantum fluctuations, e.g., charge fluctua- 
tions between half-filled t2 g and e g states in manganites 
or spin fluctuations, are not yet fully understood, while 
the origin of very low-temperature magnetism in multi- 
orbital materials remains little investigated in a realistic 
context. The hybridization-expansion continuous-time 
quantum Monte Carlo (CT-HYB) technique^— appears 
to date the most promising DMFT quantum impurity 
solver to study real materials at experimental tempera- 
tures, although most calculations so far have been lim- 
ited to high-symmetry cases or systems for which the 
hybridization function is diagonal (or almost diagonal) 
in orbital spaced " 11 ' 14 

In the present work we study the effects of com- 
monly adopted approximations on the origin of orbital 
and magnetic order in some exemplary low-symmetry 
transition-metal oxides. To do this, we use an efficient 
general implementation of the CT-HYB quantum Monte 
Carlo (QMC) LDA+DMFT solver for systems of arbi- 
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trary point symmetry and arbitrary local Coulomb in- 
teraction. In our implementation we combine a general 
Krylov 11 scheme, which we use for the low-symmetry 
cases, with the very fast segment implementation^ which 
can be used when the local Hamiltonian does not mix fla- 
vors (i.e., spin-orbital degrees of freedom). In addition, 
we use symmetric a 1 ' 1 2 to minimize the computational 
time. We present results for the orbital melting transi- 
tion in the rare-earth manganitcs i?Mn03, orbital fluc- 
tuations in the 3d 1 perovskites CaVC-3 and YTiC>3, and 
ferromagnetism in the Mott insulator YTi03. Finally, 
we investigate the regime of validity of the t 2g classical 
spin approximation often adopted to describe LaMn03 
and more general manganites. 

The paper is organized as follows. In Section II 
we briefly discuss the approach in the context of the 
LDA+DMFT method. In Section III we present ap- 
plications to rare-earth manganites, vanadates, and ti- 
tanates. We show that spin-flip and pair-hopping terms 
do not affect the Kugcl-Khomskii orbital-order transition 
and weakly affect orbital fluctuations in 3d 1 perovskites. 
We calculate the ferromagnetic transition temperature 
for the Mott insulator YTiC>3 and find excellent agree- 
ment with experiments, showing that orbital order is in- 
deed compatible with ferromagnetism in this material, 
contrarily to early hypothesis.— For LaMnC>3 we show 
that the e g two-band Hubbard model commonly used to 
study the system, in which the t 2g electrons are treated 
as disordered classical spins interacting with the e g spins 
via the Coulomb interaction, yields results in very good 
agreement with the full five-orbital 3d Hubbard model. 
Remarkably, the agreement is not only excellent for the 
occupied state in the orbitally ordered phase, but also 
very good for the orbital resolved e g spectral function 
matrix. Finally, in the Appendix we describe the details 
of our implementation of the general CT-HYB solver. 



II. MODEL AND METHOD 



Fo, F2 and F4, with £/ avg = Fq (direct Coulomb inter- 
action) and J avg = ji(F2 + F4) (exchange Coulomb in- 
teraction). In the following we find it more useful to use 
as parameters^! the diagonal element of the Coulomb 
matrix, Uo = Fq + |j, the Kanamori exchange pa- 
rameter J = I Javg and the Coulomb anisotropy 5 J = 
J{\ — + The exchange couplings for e g and 

t 2g only are then J eg = J + 36 J and J t2g = J + 5 J. We 
solve the model © with DMFT using the CT-HYB QMC 
approach as quantum impurity solver.— ~— Our implemen- 
tation of the CT-HYB QMC solver is discussed in the Ap- 
pendix. It works efficiently for systems of arbitrary space- 
group symmetry, i.e., with both a hybridization- function 
matrix and self-energy matrix in the full spin-orbital 
space. We optimize our code for modern massively paral- 
lel architectures and exploit symmetries to minimize the 
computational time. We use two approaches to calcu- 
late the trace which enters in the numerical evaluation of 
the Green function: the the segment approach^ and the 
Krylov method^ The segment approach is very fast but 
can only be used if the local Hamiltonian does not mix fla- 
vors (spin-orbital degrees of freedom). The Krylov proce- 
dure is instead general and scales linearly with the inverse 
temperature, becoming therefore particularly efficient in 
the low-temperature rcgimei 12 ' 13 Far from phase tran- 
sitions, we further enhance the efficiency by adaptively 
truncating the local trace in the Green function! 11 ' 12 Fur- 
ther details on our code are given in the Appendix. Our 
efficient implementation allows us to include in the model 
Hamiltonian ([1]) typically neglected interactions, such as 
spin-flip and pair-hopping terms or spin-orbit coupling, 
to study models with larger number of orbitals (e.g. with 
the complete 5-orbital d shell) and reach very low tem- 
peratures, as essential to study magnetic transitions. In 
the following, we use our code to systematically compare 
different models and test typically adopted approxima- 
tions on the orbital and magnetic order of a selection of 
exemplary materials. 



The most general multi-band Hubbard model for 
transition-metal oxides is given by 



H = - VV V t vl ' , ,c 1 
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Here c\ ma (c imcr ) creates (annihilates) an electron with 
spin a in orbital m on lattice site i; t mtTm , a , arc 
the hopping integrals and e mam ' a ' the elements of the 
crystal-field matrix, obtained from LDA calculations 
by constructing a localized Wannier-function basis^^ 
Umm'mm' are the screened Coulomb matrix elements, 
typically expressed in terms of the three Slater integrals 



III. RESULTS 

A. Orbital-order melting in rare-earth manganites 

The origin of the orbital-order melting transition^ in 
the rare-earth manganites i?Mn03 with the i|g e g nomi- 
nal electronic configuration has been debated since long. 
Recently^i^ we have shown that the many-body super- 
exchange interaction plays a small role in determining 
the orbital-order melting temperature Too as well as 
its trends with decreasing radius of the rare-earth ions. 
However, spin-flip and pair-hopping terms, neglected in 
previous calculations, restore the full degeneracy^— of 
the S = 1 multiplet, and could enhance the strength of 
super-exchange, or even modify the occupied orbital^* 
Furthermore previous calculations, as most many-body 
studies of rare-earth manganites, rely on the classical 
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spins approximation for ti g orbitalsj 2 ^ In such approx- 
imation the effects of the ti g spins (St 2g = 3/2) on the 
e g states is described through a local magnetic field due 
to the e,g-ti g Coulomb exchange interaction and a band- 
width rcnormahzation factor arising from the spatial dis- 
order in the orientation of the t2 g spins. However, charge 
fluctuations between tig and e g states or tig multiplct 
fluctuations, not accounted for in such a model, could 
affect the orbital-order and the occupied orbital. In this 
section we use our implementation of the CT-HYB QMC 
solver to analyze these effects. 



Role of spin-flip and pair-hopping 
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First wc analyze the role of spin-flip and pair hopping 
interactions on the orbital melting transition. The min- 
imal Hubbard Hamiltonian which is believed to retain 
the essential physics 2 ^ to study this issue is a two-band 
Hubbard model for e g states coupled to disordered ti g 
spins via the Coulomb interaction, which acts as a local 
magnetic field h = J t2 St 2 ■ Thus in Hamiltonian ([1} the 
one-electron term becomes 

t™' , , = V it U ' , 
"mam' a' u, <t,<t "mm'' 

The index m runs over the e g Wannier orbitals \x 2 — y 2 ) 
and |3z 2 — r 2 ), a z is the Pauli z matrix, while t x and 
t z are pseudospin operators acting on orbital degrees of 
freedom (t z \3z 2 - r 2 ) = l/2|3z 2 - r 2 ), t z \x 2 - y 2 ) = 
-l/2\x 2 ~ y 2 ) t x \3z 2 - r 2 ) = \x 2 - y 2 )). The ener- 
gies £jt and £t yield, repsectively, the Jahn- Teller and 
tetragonal crystal-field splitting. Finally u ayCr i = 2/3 is 
a band renormalization factor which accounts for the 
disorder in the orientations of the tig spinsi 2 ^ For the 
effective magnetic field h, wc present calculations for 
the theoretical estimate 2 ^ h ~ 1.35 eV; our results for 
the orbital-melting transition and the orbital polariza- 
tion are however weakly dependent on h in the relevant 
regime, in which e g and ti g spins arc locally aligned. For 
the e g basis, the Coulomb interaction is composed of 
density-density interactions, spin-flip and pair hopping 
terms. We use the theoretical estimates Uq = 5 eV 
and J Cg ~ J ~ 0.76 eV for the e g screened direct 
and exchange on-site Coulomb interaction! 5 ' 25 ' 26 In or- 
der to calculate the critical temperature due to super- 
exchange only, we set the crystal-field parameters to zero: 
£jt = St = 0. This disentanglement procedure has been 
proposed in Ref. 0, and was succesfully used to study 
orbital order in cuprates and manganites^i^^ 

We show in Fig. [T] the results of our calculations 
based on our CT-HYB QMC solver; we use the Krylov 
approach for the model with spin-flip and pair hop- 
ping terms and the segment method for the model with 
density-density Coulomb terms only. The figure shows 
the orbital-order transition temperature due to super- 
exchange only, T KK , for relevant elements of the series 



FIG. 1. Orbital-order transition temperature due to super- 
exchange, Tkk, versus the R 3+ ionic radius in the i?Mn03 
series (R = Dy, Tb, Nd, La). Empty symbols: Tkk (total 
energy gain) taken from Ref. Qij ; calculations were done for 
density-density Coulomb interactions and using a Hirsch-Fye 
QMC solver. Filled symbols, grey: CT-QMC (segment solver) 
and density-density Coulomb interactions only. Filled sym- 
bols, black: CT-QMC (Kryolv solver) and full Coulomb in- 
teraction. 

of rare-earth manganitcs. This figure demonstrates that 
the spin-flip and pair hopping terms affect very little the 
overall trends and even the absolute value of Tkk- These 
results all reinforce previous conclusions^ that super- 
exchange has a small influence in determining the orbital 
order to disorder transition observed in rare-earth man- 
agnitcs. 



2. Classical tig spins versus full 5-band model for LaMnOz 

Next we test the validity of the classical ti g spin ap- 
proximation for the orbital-order melting transition. To 
do this, we compare the results of the previous section 
with those obtained for the full 5-band Hubbard model 
described by Hamiltonian ([1]). To study the orbital or- 
der due to superexchange only, we again set to zero 
the crystal-field splitting within the e g doublet and ti g 
triplet; we retain however the cubic crystal field which 
splits tig and e g ; finally, we perform the LDA+DMFT 
calculations at T ~ 290 K, i.e., well below T KK . Since we 
have already shown that spin-flip and pair-hopping do 
not affect the transition temperature, we neglect them 
here to speed up the calculations. Furthermore, to com- 
pare directly the results of the two- and five-band model, 
we assume J eg -t 2g ~ h/St 2g for the e g -ti g exchange cou- 
pling and neglect other small Coulomb anisotropics. The 
LDA+DMFT calculation for the five-band model yields 
half-filled t 2g states and almost fully polarized e g states. 
The occupied e g state |0) = cos ||3z 2 -r 2 ) -sin ||x 2 -y 2 ) 
is the orbital with 8 ~ 90°, in excellent agreement with 
the results from the classical ti g spins approximation, 
which gives basically the same state. The spectral func- 
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FIG. 2. (Color on-line) Top panel: LaMn0 3 , comparison of 
the spectral function matrix obtained for the 5-band Hubbard 
model and the e g 2-band model with classical tig spins. Cal- 
culations have been performed at T ~ 290 K and Uo ~ 5 eV. 
The chemical potential fi is at energy zero for the 2-band 
model and at ~ 0.3 eV for the 5-band model. The Jahn- Teller 
and tetragonal crystal-field splittings are set to zero. Black 
line: tig spectral function. Light lines: e g spectral function 
from the 5-band model. Dark lines: e 9 spectral function from 
the 2-band model, spin up (full) and down (dashed) . The po- 
sition of the spin down Hubbard band depends on the effec- 
tive magnetic field h, i.e. on Jt 2g . Bottom panel: Comparison 
of the e g spectral function matrices, orbitally resolved. Full 
(dashed) lines: most (least) occupied orbital. Dark (light) 
lines: 2-band (5-band) model. 



tion matrix calculated for the 5- and 2-band model are 
compared in Fig. [2] This figure shows that not only 
the orbitals but also, surprisingly, the overall spectral 
function matrices are in good agreement. Because the 
five-band model includes the full dynamic of the t 2g 
electrons^ the effective Uo is larger than for the two- 
band model. By scanning different Uq between 7 eV and 
5 eV we find that Uq ~ 5.5 eV yields a gap quite close to 
that of the two-band model. This shows that in the two- 
band model the Coulomb integral Uo is screened ~ 10% 
by the t 2g electrons. The half-filled t 2g bands exhibit 
a very large gap because at half-filling the t2 g exchange 
couplings effectively enhance the effect of the Coulomb 
repulsion Uq. Finally, we find that the on-site spin-spin 
correlation function (S z a Sz a ) ~ 0.74, very close to the 
value 0.75 expected for aligned e g and S t2g = 3/2 t 2g 
spins. For what concerns the sign problem, we find it 
negligible for all of these calculations (the average sign is 
~ 0.99 in the worse case). 
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n 2 


nz 


CaV0 3 


0.47 


0.28 


0.25 


YTiOs 


0.98 


0.01 


0.01 



TABLE I. Occupations m of the natural orbitals (with m > 
fii+i) at T = 190 K in CaV0 3 and YTi0 3 obtained by diag- 
onalizing the occupation matrix. 



Orbital fluctuations and magnetism in CaV0 3 
and YTi0 3 



The importance of orbital fluctuations in the physics of 
3d 1 perovskites has been debated since long£ i 15 i 16 i 27 ~ — 
Single-site DMFT calculations have shown that in the 
presence of crystal-field splitting Coulomb repulsion 
strongly suppresses orbital fluctuations.— However, these 
conclusions were based on a Hubbard model with density- 
density Coulomb interactions only. In this section we an- 
alyze the effect of the neglected spin-flip and pair-hopping 
Coulomb interactions. Furthermore, exploiting our effi- 
cient CT-HYB solver, we address the issue of the nature 
of the low temperature (30 K) 15 ' 30 ferromagnetic transi- 
tion in YTi0 3 . 



1. Orbital fluctuations 

The minimal model to consider for 3d 1 transition-metal 
oxides is a three-band Hubbard model for the t 2g bands 
including spin-flip and pair hopping terms, and with 



-mom' a' 



t 



mam' a' 



t U ft , 
b m.m' u a,(J' 



where m,m! = xy,xz,yz. For the Coulomb parame- 
ters we use Uq = 5 eV and J t2g ~ 0.68 eV (CaVOa) 
or J t2g = 0.64 eV (YTi0 3 ) from theoretical estimates 
and previous works^i Because the local Hamiltonian 
mixes flavors even in the crystal-field basis, we perform 
the LDA+DMFT calculations using the Krylov version 
of our general CT-HYB QMC solver. 

In Table U we show the occupations rij of the natural 
orbitals at -190 K in CaV0 3 and YTi0 3 . We find that 
CaV0 3 is a paramagnetic metal with a small orbital po- 
larization. Instead, YTi0 3 is a paramagnetic insulator 
with orbital polarization p = n\ — (n- 2 + ?i 3 )/2 — 1, i.e. 
basically full (orbitally ordered state). For this system, 
the double occupancies at 290 K are small, i.e., we find 
lEm^v^v) ~ 0.015 for YTi0 3 . The occu- 
pied orbital is 0.611|xy)-0.056|a;z}+0.789|j/z). We find 
occupied state and orbital polarization are basically the 
same with full Coulomb and density-density approxima- 
tion. Previous calculations^ in which spin-flip and pair- 
hopping terms have been neglected and T — 770 K are 
in line with these results. This shows that spin- flip and 
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pair-hopping terms do not change the conclusion that 
orbital fluctuations are strongly suppressed in the Mott 
insulator YTi0 3 . In the CT-HYB QMC simulations the 
average sign is ~ 0.9 for YTiC>3 and ~ 0.95 for CaV03. 



2. Ferromagnetism in YTiOz 

YTiC>3 is one of the few ferromagnetic Mott insula- 
tors. Neutron scattering experiments pointed out early 
on the difficulties in reconciling ferromagnetism and the 
expected orbital order,— and there have been sugges- 
tion that the ferromagnetic state could rather be asso- 
ciated with a quadrupolar order and large scale orbital 
fluctuations.— However, second-order perturbation the- 
ory calculations indicate that ferromagnetism and orbital 
order could be reconciled, provided that the real crystal- 
structure of YTiC>3, including the GdFe03-type distor- 
tion (tilting and rotation of the octahedra, and deforma- 
tion of the cation cage) is taken into account^ To clar- 
ify this point, we check the instability towards ferromag- 
netism of the three-band tig Hubbard model obtained 
for the experimental structure of YTiC>3. With this ap- 
proach we calculate the ferromagnetic transition temper- 
ature Tc due to super-exchange alone in the orbitally 
ordered phase. Since experimentally Tc ~ 30 K, we have 
to perform LDA+DMFT calculations down to very low 
temperatures, which becomes possible with the CT-HYB 
QMC solver. On lowering the temperature, we find that 
the sign problem becomes sizable (average sign ~ 0.7 
at 40 K). However, we can basically eliminate it (average 
sign ~ 0.97) by performing the LDA+DMFT calculations 
in the basis which diagonalize the crystal-field matrix, 
even though the hybridization function has off-diagonal 
terms of comparable size in the two bases. In Fig. [3] we 
show the LDA+DMFT magnetization m(T) of the t 2g 
states as a function of the temperature. Remarkably, we 
find a transition at about 50 K, in excellent agreement 
with experiments^ which yield Tc ~ 30 K; the overes- 
timation can be ascribed to the mean-field approxima- 
tion, and to the fact that, since the critical temperature 
is very small, it is sensitive to tiny details. The occu- 
pied orbital does not change significantly in the magnetic 
phase, indicating that the occupied orbital remains the 
one that diagonalizcs the crystal-field matrix, i.e., in the 
magnetic phase there is no sizable change of orbita&i2 
due to super-exchange. 




T/K 

FIG. 3. Ferromagnetic spin polarization as a function of tem- 
perature in YTi03. The plot shows a transition at the critical 
temperature Tc ~ 50 K, slightly overestimating the experi- 
mental value Tc ~ 30 K, as one might expect by mean-field 
calculations. 



(i.e., if the local Hamiltonian does not mix flavors) we 
use the fast segment solver. In more realistic situations 
we use the Krylov approach and, away from phase tran- 
sition, trace truncation. We find that in all considered 
cases the minus sign problem mostly appears when off- 
diagonal crystal-field terms are present (and is strongly 
suppressed in the basis of crystal-field states), while off- 
diagonal terms of the hybridization function matrix are 
not as critical^ We show that spin-flip and pair-hopping 
terms hardly affect the strength of the super-exchange 
orbital-order transition temperature in rare-earth man- 
ganites. We show that the classical tig spin approxima- 
tion for LaMn03 works excellently, not only for what 
concerns orbital order, but, surprisingly, also for the 
overall shape of the spectral function matrix. We show 
that spin-flip and pair-hopping terms also do not change 
the conclusion that orbital-fluctuations are strongly sup- 
pressed in YTi03. Furthermore, we calculate the critical 
temperature for ferromagnetism in the orbitally ordered 
phase, and find excellent agreement with experiments. 
This shows that that the predicted orbital order is fully 
compatible with ferromagnetism. 



IV. CONCLUSIONS 



We implement an efficient general version of the 
continuous-time hybridization expansion (CT-HYB) 
quantum Monte Carlo solver, which allows us to in- 
vestigate ordering phenomena in strongly correlated 
transition-metal oxides in a more realistic setting. Our 
implementation of CT-HYB QMC works for systems of 
arbitrary symmetry In cases where symmetry allows it 
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Appendix A: General CT-HYB solver 

In this appendix we fix the notation and explain 
the details of our implementation of the general CT- 
HYB quantum-impurity solver. The DMFT quantum- 
impurity Hamiltonian is H = H\ oc + i?bth + -ffhyb, where 



qq act 



-Hbth = E e 7 & 7 6 7 , 

7 

#hyb = 51 51 [^7,« c a b 7 + ^ c -] 



The combined index a = ma labels spin and orbital 
degrees of freedom (flavors). For the bath we use, 
without loss of generality^ the basis which diagonal- 
izes -ffbth, with quantum numbers 7. Finally, we define 
£aa = Ea&~ Aejj, where £ a& is the crystal- field matrix 
and Ae^g is the double counting correction; in the cases 
considered in the present paper the latter is a shift of the 
chemical potential fi. 



where c^\t) = e T ( ff i--^) c (t) e -T(ff loc - M jV) and N is thc 
total number of electrons on the impurity. For expan- 
sion order m = 2n, the vector a = (a\, a.^ ■ . ■ a n ) gives 
the flavors on associated with the n annihilation oper- 
ators on the impurity at imaginary times Tj, while the 
a = (0:1, a.2 . ■ ■ a n ) are associated with the n creation 
operators at f,. The second factor is the trace over the 
non-interacting bath, which is given by the determinant 



4"L(t,t) = det[i^, (t,t)] 
of the n x n square hybridization-function matrix with 



(n) 



matrix elements [-^a a( T ) T )]« 
by 



FoutoLiin' - n) given 



1 + e~P e -< 

7 



-e"^ T t > 

e -M/3+r) T < Q _ 



On the Fermionic Matsubara frequencies, w n , its Fourier 
transform 



1. Hybridization-function expansion 

By expanding the partition function in powers of 
i?hyb and going to the interaction picture _ffhyb('?") = 
e r(H btb +H }oc ) Hhyhe -T(H bti ,+H loc ) with p = i/k B T we ob- 
tain the series 



Tr \eT^ Hb ^ +H ^Te~ $ drH hyb{r)\ 

OO „ 



M 

dr Tr T 



-/3(ffbth+i/loc 



I #hyb(l 



where T is the time order operator, r = (ri, T2, . . . r m ) 
with r, + i > Tj and 



(m) 



In the trace only terms containing an equal number of 
creation and annihilation operators in both the bath 
and impurity sector, i.e., only even expansion orders 
m = 2n contribute. Introducing the bath partition func- 
tion Zbth = Tr e~^ ffbth , the partition function can be 
factorized 



00 /•(«) /•(«) . , 
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with 



The first factor is the trace over the impurity states 
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-f3(H loc -nN) 



is related to the bath Green-function matrix Q by 

Faction) = i^n^aa — ^aa~ (G)aa( 0J n), 

as can be shown by downfolding2& 
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to the impurity block (i = 0). Here the matrix elements 
of H and I are given by (H ) a a = 4a and (Jo)aa = 
S at& , while (V 0t i) m = V &ti , and (V ifi ) ia = V itCx . 

To speed up the calculations, we exploit symmetries. 
If Nt blocks of flavors are decoupled by symmetries, the 
hybridization function matrix is block-diagonal in those 
flavors. We then write the partition function in terms of 
the expansion orders n& in each block, with n = X^ti n bt 



Zbth 

with 



N b 00 

nE 

,6=1 n b =0 



(n b ) 

dn 



(n b ) 

df b 



E 

otb<x b 



An) 



k(r,r) 



4:L(T,r) = n4;l(n,n) 



and 



i n) a (^^) = Trr 



b=l 



-P(H loc -nN) 



N b 



6=1 i—rih 
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FIG. 4. Convergence of the Krylov approximation \ip(T)) r 
to \t^(r)) = e~ ( - Hl °"~ E °^ T \il)) for a representative test case (5- 
orbital model, half filling). The figure shows the difference 
A(r) = | \4>{t)) t — I'M 7 ")) I- Symbols (in order of increasing 
size): r = 0.005, 0.05, 0.5, 5 and 100. 



2. Segment solver and Krylov approach 

Calculating the trace over the impurity states involves 
propagating states in the impurity Hilbert space. For 
models with many orbitals this can become very demand- 
ing. We therefore use a multi-approach scheme. When 
the on-site Hamiltonian conserves the flavors we use the 
so-called segment approach,— which is extremely fast. In 
such cases only terms with an equal number of creation 
and annhilation operators per flavor contribute to the lo- 
cal trace, and it is convenient to express the partition 
function in expansion orders n a for flavors a. The parti- 
tion function then can be rewritten as 



^bth 



IIE 

o=l n a =0 



("a) 



a)_ 

dr a 



» 



k(r,r). 



Here r = Yla=i Ta an< ^ T = J2a=i r ai while the vectors 
a = ^2 a a and a = Y] ° 1 a a have the n a compo- 
nents a„i = ctai = a. The local trace factors into 



Na 1 

a—1 i—n a 




,-12 aa >[&aa-n)Sa,*' + k**a>]la 



where l aa i is the length of the overlap of the r segments 
a and a', s a = sgn(r a i — f a i) is the Fcrmionic sign, and 

Uaa* = Uaa'a'a + V aa' aa' IS the interaction. 

In all the cases in which the local Hamiltonian mixes 
flavors, we adopt the Krylov method*^ At the begin- 
ning of the DMFT loop we calculate all the eigen- 
states of H\ oc . and their energies {E n }; a given 

state \*f? n ) is then propagated with e~ TlEn ; the first 
creation or annihilation operator met generates a new 



state |\E'), which we propagate with e -( T 2- T i) ff ioc obtain- 
ing | x I'(t2 — Ti)); we repeat the procedure till the last 
creation or annihilation operator is met. At the core 
of the procedure are the matrix-vector multiplications 
and the propagation of vectors. For the first aspect, 
we work in the occupation number basis, in which H\ oc , 
and the creation and annihilation operators are sparse 
matrices. Additionally, we arrange the states according 
to the symmetrie s 10 ' 12 of H\ oc , so that we have sparse 
block-diagonal matrices and can exploit to the maxi- 
mum efficient sparse-matrix multiplication algorithms. 
We find that this typically reduces the CPU time by, 
e.g., about 15% for a three-band model. We use the 
Krylov approach to calculate \^(t)) = c~ HlocT \ty). First 
we construct the Krylov space of order r, JC r (\^)), i.e., 
the space spanned by |#), H loc \^), H^JV) . . . #£J*). 
By means of the Lanczos^ technique we construct 
an orthonormal basis for /C r (|\&)), {|fc)}; in this basis 
ifioc is tridiagonal with eigenstates {\l}} and energies 
{si}. The matrix exponential e~ HlocT is approximated 
by its projection onto the Krylov space, e~ HlocT \' l b) ~ 
|\&(r)) r = X)[=o e_Tei IO(^l v l')- This procedure converges 
very rapidly with r, typically for r much smaller than 
the dimension of the Hilbert space ; 37 ' 38 as illustrated in 
Fig.SJ We find that the convergence slightly deteriorates 
with increasing r and the complexity of the Hamiltonian 
(realistic Coulomb vertex, crystal- field matrix) , but typi- 
cally 2-3 steps are sufficient to obtain accurate results. To 
best exploit the power of the method, we keep r flexible. 
Furthermore, to avoid that the norm of the state becomes 
very large during the propagation, we set Eq to zero, 
i.e., substitute e~ rHloc with e ~ T ( Hloc ~ Ea \ In addition, 
the procedure (propagation and creation/ annihilation) is 
carried out from both the left and the right side of the 
trace, to minimize the work needed to measure, e.g., the 
Green function matrix. Finally, at low temperatures or 
far from phase transitions we use the eigenvalues of H\ oc 
to determine the relevant energy window and truncate 
adaptivcly the outer bracket of the trace. This further 
reduces the CPU time. 

The performance of our CT-HYB QMC solver (Krylov 
and segment version) on the Julich BlueGene/Q, and 
comparison with Hirsch-Fye QMC, is shown in Fig. [5] 



3. Green-function and occupation matrix 

The partition function (|A1[) can be seen as the sum 
over all configurations c = {a,Ti, otifi, n} in imaginary 
time and flavors. In a compact form 

z = E (Z)c = E w - - E si e n ( Wc ) . 
c c :• : 

where in the last term the sum is over a sequence of 
configurations {c} sampled by Monte Carlo using |u> c | as 
the probability of configuration c. In the segment solver 
approach, we parametrize the configurations by intervals 
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FIG. 5. (Color on-line) Scaling of our CT-HYB QMC 
LDA+DMFT code on BlueGene/Q. Black line: Hirsch-Fye 
Quantum Monte Carlo solver, 2 orbitals. Other lines: CT- 
HYB Krylov (dark) and CT-HYB segment (light). Symbols: 
Two- (circles), three- (triangles) and five- (pentagons) band 
model. Open symbols: Truncated local trace. All points cor- 
respond to calculations of high quality (and with comparable 
error bars) for the systems considered in this work. 



[0,/3) (time-line), occupied by a sequence of creators and 
annihilators, which define segments on the time-line. The 
basic Monte Carlo updates are addition and removal of 
segments, antisegments or full lines£ In the Krylov solver 
approach we use the insertion and removal of pairs of 
creation and annihilation operator o 9 ' 10 as basic updates. 
In addition, we shift operators in time^^ and exchange 
the configurations of blocks or flavors^, (global moves). 
Finally, a generic observable O can then be obtained as 
Monte Carlo average 



O 



E{ C } si gn(> c ) 



where (0) c is the value of the observable for configura- 
tion c, and c runs over the configurations visited with 
probability \w c \ during the sampling. The average ex- 
pansion order increases linearly with the inverse temper- 



ature. For the case of YTi03, at ~ 40 K, the average 
expansion order is n ~ 40. 

We calculate the Green function matrix in two ways, 
directlj&iS and via Legendre polynomials^ In the first 
approach, the Green function matrix is obtained as 
Monte Carlo average with (0) c = (G a a)c, and 

N b n b 

{G a& )c = Y^ E ^{ T ^n 3 -f bl )[M {nb) ] bjM 8 ab]a 8 ab%a . 

6=1 f,3=l 

Here = [F^]" 1 is the inverse of the hybridization- 

function matrix, which we update at each accepted move, 
while A is given by 

A(r,rO = -i| 6iT ~ T,) T ' >0 
P\-S{t-(t j +0)) r'<0 

and the <5-function is discretized. In the second approach 

we c 

with 



we calculate the Legendre coefficients (0) c = (G l aa ) c , 



N b n b 

EE 

6=1 i,j=l 



(G l aa )c = E E P ^ ~ r M )[M^\ 3M 5 abDa 5 &b ^ 



Pi{r) 



V2TTT 



Pl(x(r)), t>0 
-Pi(x(t + 0)), r<0 



where pi{x) is a Legendre polynomial of rank I, with 
x(t) = 2t//3 — 1, and we reconstruct the Green func- 
tion matrix from 



1=0 



V21 + 1 



Pi(x(T))G l a£ 



For what concerns occupations, in the segment solver we 
calculate them from the total length of the segments of 
the different flavors^ in the Krylov solver we obtain them 
in two ways, directly from the Green's function and by 
explicitly inserting the occupation number operator at 
the center of the operator sequence (r = /3/2)and calcu- 
lating the corresponding trace.^i The off-diagonal ele- 
ments of the local occupation matrix (cJ^Cg), which can- 
not be obtained by inserting the corresponding operators 
at t = (3/2,— are extracted from the Green function ma- 
trix only. 
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